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Abstract 

Genotype-phenotype (GP) maps specify how the random mutations that change genotypes generate variation by altering 
phenotypes, which, in turn, can trigger selection. Many GP maps share the following general properties: 1) The total number 
of genotypes N G is much larger than the number of selectable phenotypes; 2) Neutral exploration changes the variation 
that is accessible to the population; 3) The distribution of phenotype frequencies F l ,=N p /N G , with N p the number of 
genotypes mapping onto phenotype p, is highly biased: the majority of genotypes map to only a small minority of the 
phenotypes. Here we explore how these properties affect the evolutionary dynamics of haploid Wright-Fisher models that 
are coupled to a random GP map or to a more complex RNA sequence to secondary structure map. For both maps the 
probability of a mutation leading to a phenotype p scales to first order as F p , although for the RNA map there are further 
correlations as well. By using mean-field theory, supported by computer simulations, we show that the discovery time T p of 
a phenotype p similarly scales to first order as l/F p for a wide range of population sizes and mutation rates in both the 
monomorphic and polymorphic regimes. These differences in the rate at which variation arises can vary over many orders of 
magnitude. Phenotypic variation with a larger F p is therefore be much more likely to arise than variation with a small F p . We 
show, using the RNA model, that frequent phenotypes (with larger F p ) can fix in a population even when alternative, but 
less frequent, phenotypes with much higher fitness are potentially accessible. In other words, if the fittest never 'arrive' on 
the timescales of evolutionary change, then they can't fix. We call this highly non-ergodic effect the 'arrival of the frequent'. 
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Introduction 

Darwin's account of biological evolution [1] stressed the 
importance of natural selection: If some individuals are better 
adapted to their environment than their competitors, their 
offspring will come to dominate the population. The fittest survive 
and the less fit go extinct. Yet selection alone is not sufficient to 
drive evolution because natural selection reduces the very 
variation that it requires to operate. It was only recognised well 
after Darwin's day [2], in part through the success of the Modern 
Synthesis, that the fuel for selection is provided by mutations that 
make offspring genetically different from their parents. Crucially, 
mutations change genetically stored information (the genotype) while 
selection operates on the physical expression of this information 
(the phenotype). Understanding the relation between genotypes and 
phenotypes - the GP map - is therefore crucial to understanding 
evolutionary dynamics [3]. 

GP mappings have been studied at different levels of abstraction 
[4] The most basic systems are concerned with the sequence-to- 
structure(-to-function) relation of single molecules such as RNA [5] 
or proteins [6-8], but higher-level systems such as protein 
complexes [9], gene-regulatory networks [10] and developmental 
networks [1 1] have also been studied. Even though these GP maps 



arise in quite different contexts, they share several interesting 
properties: 

1) Most basically, the number of possible genotypes Ng is 
typically much greater than the number of possible pheno- 
types Np, so the map is many-to-one. As a consequence, many 
mutations may conserve the phenotype, leading to mutational 
robustness. Important prior work has linked such robustness 
to the concept of neutral spaces, namely the set of all 
genotypes that map to a particular phenotype, with the 
additional property that they be linked by neutral mutations 
[4,5,12]. 

2) Even though Np«Ng, the accessible genetic neighbourhood 
of a single genotype g that generates a given phenotype p may 
include significantly fewer alternative phenotypes (potential 
variation) than is found in the neighbourhood of the (neutral) 
set J\f p of all Np = [N~ p \ genotypes that map onto phenotype p. 
Exploration of a neutral space can therefore increase the 
variety of phenotypes discovered by a population [13,14]. 

3) Perhaps the most striking commonality of these GP maps is a 
strong bias in assignment of genotypes to phenotypes: Most 
phenotypes are realised by a tiny proportion of all genotypes, 
while most genotypes map into a small fraction of all 
phenotypes. This property is shared by all the GP maps we 
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noted before. Typically the number N p of genotypes per 
phenotype p and the related phenotype frequencies 
F p =N p /Ng can vary over many orders of magnitude. Such 
huge variations are likely to have an effect on the course of 
evolution. 

In this paper we study the evolutionary dynamics of a classical 
Wright-Fisher model, but with explicit microscopic GP maps that 
capture the three generic properties of such maps introduced 
above. Motivated by the strong bias in the distribution of the F p 
observed for many GP maps, we derive a mean-field like 
approximation for the average probability § pq that a mutation 
will change a genotype that generates phenotype q into one that 
generates phenotype p. This approximation greatly simplifies the 
dynamics, allowing us to calculate analytic expressions for 
quantities such as the median time T p for phenotype p to first 
appear in the population as a function of population size N, the 
point mutation rate p., genome length L and the mutation 
probabilities <f> pq . 

These approximations are then tested against extensive 
simulations of two models: firstiy, a simple GP map where the 
genotypes are randomly assigned to phenotypes according to a 
pre-determined distribution for the frequencies F p and secondly, 
the well-known mapping of RNA sequence to secondary structure 
[4,5,15], which is more complex, but also more biologically 
realistic. We focus on the case where a population of N individuals 
has initially equilibrated at a fitness maximum given by phenotype 
q, and then measure the median time T p for alternative phenotype 
p to first arise in the population. 

Our analytic expressions agree quantitatively with the simula- 
tions in the polymorphic limit where NLp » 1 , and also in the 
opposite monomorphic limit NLp « 1 . In between these regimes a 
single scaling factor must be included. In all regimes the median 
discovery time T p cc\/<j> pq . For the random model (j> pq ~F p ; this 
scaling also holds for the more complex RNA mapping, although 
there is significantly more scatter due to local correlations within 
the neutral spaces and for some phenotypes we find (j> pq = 0 even 
though F p is large (this can be due to biophysical constraints 
explained for example in ref. [16]). Despite such higher order 
effects, the variation of the F p over many orders translates directly 
into the T p . More frequent (higher F p ) phenotypes are therefore 
discovered more rapidly and more often along evolutionary 
trajectories. In this way the structure of the GP map can play a key 
role in determining evolutionary outcomes. 

Finally, we employ the RNA GP map to study the case where 
two phenotypes pi and p2 are both more fit than the source 
phenotype q, but where F p \ »F p 2 (or more accurately <l> p i q »<l ) p 2q)- 
Direct simulations show that phenotype pi, which is more 
frequent, is much more likely to fix in the population, even if its 
fitness is much lower than that of p2, an effect we call 'the arrival 
of the frequent'. 

Results 

Theoretical framework 

We study the evolution of a population of N asexual haploid 
individuals. Each individual i carries a genotype gi of L letters 
taken from an alphabet of size K. The individual's phenotype pi is 
determined from g, via the GP map. The population evolves in in 
discrete, non-overlapping generations according to the classical 
Wright-Fisher model for haploid individuals: At each generation 
T, N parents are drawn with replacement with probability 
proportional to their fitness 1 + with the constraint that the 



population size (or carrying capacity) N is fixed. Each parent gives 
rise to one offspring, and the offspring make up the population for 
the next generation. During reproduction, each base in the 
genotype of length L mutates to a random alternative base with 
probability p. The number of mutations (that is, the Hamming 
distance) d between parent and offspring is thus distributed 

binomially according to h(d) = (^^jp. d {l — p) L ~ d . In this way the 

set {gi} of N genotypes changes at each generation. 

The expected number of individuals with phenotype p that 
arises at generation t can be written as: 

N L 

m p (t) = Yl h(d)%(gi,s„d) (1) 
i d=\ 

where <S> p (gi,Si,d) is the probability that a d fold mutation of 
genotype gi (selected for reproduction according to fitness 1 + s,-) 
generates an individual with phenotype p. It takes into account the 
mutational connections between the Nq = K L genotypes that 
make up the GP map. The probability of not finding p is 
approximately given by the Poisson distribution as exp( — m p (f)). 

While exact, these dynamic expressions depend implicitly on 
time through stochastic changes in the set {gi}, and are typically 
very hard to solve. In order to gain intuitive insight, we employ a 
number of simplifications and approximations, motivated in part 
by the general properties of GP maps discussed in the 
introduction. First, we assume that Lp«l so that for d>l, 
h(d)«h(l)»Lp, which means that we can ignore higher order 
mutations (terms with d>l in Eq. (1)). For a given source 
phenotype q (where the fitnesses of all genotypes mapping into q 
are equal, and so we take 1 + S q = 1 for simplicity) we can then 
calculate the mean probability (j> pq that a single point mutation will 
generate another phenotype p: 

where the sum is over the set N q of all N q genotypes that generate 
phenotype q (see also Figure 1). It is convenient to introduce the 
robustness of phenotype q as the average probability over all M q of 
neutral mutations: p = <j> qq - If we consider the case where at 
generation t—l the whole population is on N q , then Eq. (1) 
simplifies in this mean-field (or pre-averaged) approximation to: 

m p {t) = NLpi/ pq (3) 



The polymorphic limit. It NLp » 1 then the population naturally 
spreads over different genotypes, a regime called the polymorphic limit. 
Consider the case where l+s p =S qp so that the population remains on 
A/^which is one way to model neutral exploration. In the mean-field 
approximation the expected number of individuals with phenotype 
p produced per generation is now independent of time, and given 
by Eq.(3), as long as double mutations can be ignored. The time 
T p {u) when on average the probability of having discovered p is a 
(so that the median discovery time of/? is T p (l/2)) is then given by: 
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Figure 1. Illustration of the mean field approximation. A) An example genotype space: Each point corresponds to a unique genotype; shape 
and color of the marker indicate the phenotype. Genotypes joined by edges can be interconverted by single mutations. Edges for neutral mutations 
share the color of the (conserved) phenotype, non-neutral mutations are shown as black dashed lines. The shading of the genotypes illustrates the 
number of individuals carrying the respective genotype in a hypothetical population. The mutations away from the genotypes occupied by the 
population determine the accessible phenotypes. 8) Our meanfield approximation averages over the internal structure of neutral spaces. So neutral 
spaces are represented by the markers of their phenotypes only, with the size representing the neutral space size (ie. number of genotypes in the 
space). The uniform shading of the blue neutral space implies that in the meanfield approximation, the population is assumed to continually explore 
the neighbourhood of its entire neutral space. Mutational outcomes are thus determined from the local frequencies of phenotypes around the 
neutral space, as measured by the <j> M coefficients. This mean field approximation allows us to derive analytic forms that can be compared to 
simulations of the full GP map. 
doi:1 0.1 371 /journal.pone.0086635.g001 



Tp(0£) = 



-log(l-o() 
NLf4 



(4) 



Eqns. (3-4) should provide a good approximation of the full 
dynamics in the limit that N is large enough that variations 
between individual genotypes gveA/ q are averaged out, in other 
words, for the case where the 1 -mutant neighbourhood of the 
population is similar to that of the whole neutral space. 

The monomorphic limit. Neutral spaces can be astronomically la- 
rge [1 7], much bigger than even the largest viral or bacterial populat- 
ions. In that case, the local neighborhood of the population may not 
be fully representative of the neighborhood of the entire space. This 
scenario can most easily be understood in the monomorphic limit 
where mutants are rare, NLp«l, and exploration is dominated 
by genetic drift. Every neutral mutant has a probability of 1 /N to 
go to fixation, allowing the population to move to a new genotype. 
Thus the timescale of fixations is Kimura's famous result [18] 
Zf = \j(Lfip), where the robustness p is the probability that a 
mutation is neutral, so that Ljip is the rate of neutral mutations. 

Between fixations, the population undergoes periods of geno- 
typic stasis in which only the 1 -mutant neighborhood of the 
current genotype g is explored by (rare) mutations. As there are 
(K—l)L adjacent genotypes, the timescale of this exploration is 
x e = (K- \)L/{NLfi) = (K-l)/(Nfi). 

It is instructive to compare the ratio £ of these two time-scales, 
defined via 



Tf N 

x~ e ~ (K-l)Lp * L 



(5) 



We can use this dimensionless ratio to distinguish between 
different dynamic regimes. If C»l, fixation takes much longer 
than exploration. If we define n% as the number of local 
neighbours of the genotype g mapping to phenotype p for the 



current population, then in this limit, phenotypes with n g p > 0 are 
produced continuously (on a time-scale given by x e ) until the 
population moves to a different genotype. The dynamics under 
strong genetic drift therefore induce short-term correlations in the 
mutant phenotypes. Since C,~N jL, we call this regime the large 
population limit. 

In the opposite extreme £ « 1 , which we call the large genome 
limit, the population typically moves to a different genotype before 
all accessible mutants have been explored. In this regime, we do 
not expect short-term correlations in the mutant phenotypes, 
simply because every mutant occurs only very rarely. 

Actual discovery and neutral fixation times can show strong 
fluctuations. As our evolutionary process is a Markov process - the 
next set of mutants depends only on the parents, not on earlier 
mutants - the first discovery time of a neighbour genotype as well 
as the arrival time of the neutral mutant "destined" to be fixed, are 
distributed geometrically (or exponentially in a model with 
continuous time). Thus the mean of x e or Xf is equal to the 
respective standard deviation, and any particular evolutionary 
trajectory can be very different from the average behaviour. 

Let i be the actual time the population stays at the current 
genotype. In the continuous time approximation, x is distributed 
exponentially with mean Ty. If the genotype g has n g p mutations 
leading to p then the probability that p is found during this time is 
1 — exp( — n g p x/x e ). Integrating over the distribution of x, we have 
the probability P(j&) that phenotype p is discovered before the 
next neutral fixation: 



,co dx . 
o V 



-(1- 



> _t/ V = l. 



1 



i+«K 



(6) 



If fixations are the rate-limiting step (ie. £ » 1), P— > 1 if n g ^ 0, as 
each neighborhood is searched exhaustively before the population 
moves on. On the other hand, if fixation is faster than exploration 
(£« 1), the introduction of alternative phenotypes is determined by 
random fluctuations, as most available mutants are not produced. 
To leading order, we find P(npxn^ = Nnf/((K-l)Lp). We 
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note that the inverse dependence on p arises from Xf. More robust 
neutral spaces are explored faster, but therefore less thoroughly. 

The dynamics in the monomorphic regime are thus relatively 
straightforward. But whether some new phenotype p is discovered 
still depends on the structure of the neutral space which in turn 
determines how the available phenotypes change upon a neutral 
fixation. To describe this structure, we turn again to a mean-field 
approximation: The mutational neighborhood of each particular 
genotype geJ\f q resembles the average over Af q . As the mean 
number of mutations per genotype leading to p is given by 
n pq = {K—X)L(j> pq , the probability that p is accessible after a 
neutral fixation is 1 — exp( — h pq ) x n pq (the approximation is valid 
provided n pq « 1 , that is p is not accessible from every genotype in 
the source neutral space; of course, this is just the condition we are 
interested in, as otherwise neutral exploration would not typically 
be necessary for phenotype p to arise). 

Over a large number of generations (t»t/), a monomorphic 
population explores its neutral space uniformly [19]. Assuming 
that n s > 1 can be ignored in practice, we have 
T p (a) = —Xf\og(\—v)/{n pq P(Y)). The first discovery time in the 
large population limit becomes: 



7» = 



-T/log(l-a) 



-log(l-a) 



n pq L 2 (K-l)pp<j> pq 
whereas in the large genome limit we obtain 



t p (*)- 



-x e log(l-a) 



-log(l-a) 

NLp4 pq 



(7) 



(8) 



which has the same form as the polymorphic limit, Eq. (4): When 
the population is too small (compared to the genome length), the 
exploration of each genotype's mutational neighborhood is 
typically incomplete. Then, just as in the polymorphic limit, only 
random fluctuations determine which accessible genotypes are 
actually realized by the population. 

Finally, let us compare our results for large populations in the 
monomorphic and polymorphic limits. Most importantly, in both 
cases T p is inversely proportional to <j> pq : Rare phenotypes are hard 
to find. Comparing Equations (4) and (7), the only difference is 
that N in the polymorphic regime is replaced by L(K—\)p in the 
monomorphic limit. This difference is intuitive: When the 
population is diverse, every new individual helps exploration and 
reduces discovery times. But if all individuals have the same 
genotype, simply having "more of the same" does not make 
neutral exploration faster. However, repeated mutants may 
influence the fixation of adaptive phenotypes. 

These results suggest that for intermediate NLp there should be 
a smooth transition between these two regimes. To quantify the 
crossover we introduce a factor y that multiplies N in Eq.(4); we 
expect that y-»l as either NLp becomes very large (the 
polymorphic limit) or N«L (the large genome limit), and that 
y^(K—\)Lp/N as NLp«\ and N»L (the large population 
monomorphic limit). 

Simulations in model GP maps 

In order to test our mean-field theory we study two kinds of GP 
maps that both include the generic properties of GP maps that we 
introduced earlier. 

Random GP map. In the random GP map, the total number 
of phenotypes Np and the frequencies {F p } can be set arbitrarily 
(subject to the normalization constraint ^2pL\F p = l). The 



K L x F p genotypes mapping into phenotype p are distributed 
randomly in genotype space. The statistical properties of the map 
are thus determined by the parameters L, K, and the set {F p }. 

Studying this map has two motivations: First, ignoring some 
biophysical detail may help illuminate generic features shared by 
the systems described in the introduction. Second, a simple model 
may clarify which deviations from our theory arise from 
population dynamic effects rather than from detailed (and 
system-specific) structure in the GP map. 

In this simple model, correlations between genotypes are absent, 
facilitating analysis of the resulting neutral spaces. For example, 
<j) pq =F p is a good approximation as long as Np « Nq and 
N q ,N p »l. Also, there is a percolation threshold 
X(K)=l -K- ] ' (K - ]) : thus only phenotypes with F q > l(K) have 
completely connected neutral spaces [20]. 

Here we study a particular random GP map with L= 12, and 
K = 4 (as in DNA and RNA) so that there are 
N G = 4 l2 x 1.68 x 10 7 genotypes. These map onto jV/> = 58 
phenotypes distributed with frequencies F p cc\.2~ p . The F p vary 
over about 5 orders of magnitude, a range similar to the F p of 
L = 12 RNA (see also Figure SI). To make sure that the largest 
neutral space percolates, its frequency is set separately as 
F\ = 0.5 >/t(4) = 0.37. For several values of p., we simulated 
iV=1000 individuals for up to 7 x 10 10 generations. The fitness 
was set as 1 +s p = 5 Pi \ so that we are effectively modelling neutral 
exploration on the space A/i , which is convenient for measuring all 
T p . We measured first discovery times for the 57 alternative 
phenotypes over 100 independent simulations to obtain the 
median time T p . 

Figure 2A depicts these median discovery times T p for 
simulations ranging from the polymorphic regime NLp » 1 to 
the monomorphic limit NLp« 1 (see also Figure S2). We note the 
following: 

1) For all regimes the T p vary over many orders of magnitude, 
but they are found in fewer generations for larger p. 

2) Locally frequent phenotypes (i.e. those with high <j) pq ) are 
much easier to discover. The inset of Figure 2A shows that 
<j> pq »F p , so this conclusion carries over to frequent pheno- 
types with large F p . 

3) A subset of the phenotypes with <j) pq ><j) L =l/(K—l) 
Lx 0.028 are likely to be in the one-mutation neighbourhood 
of any genotype. In the monomorphic regime these are are 
then found by exploration of a genome so that T p is given by 
Eq. (8), which has the same form as the polymorphic limit, Eq. 
(4), as can be seen in Figure 2A. Discovery times cross over to 
the regime where neutral exploration is required when 
(j> pq «(l>i. Such behaviour can be viewed as a finite size effect: 
Np typically increases with L. Therefore the largest F p will 
likely decrease for larger systems, so that a smaller fraction of 
phenotypes can be found without neutral exploration. 

4) In the fully polymorphic regime where each individual 
essentially explores independently, any phenotype with 
(j> pq > 1 /(NLp) is likely to be part of immediately accessible 
standing variation [21] in the initial population, and is therefore 
found quickly. Indeed, in Figure 2A for p=\0~ 2 , where 
NLp = 120, these phenotypes are typically found in one or 
two generations on average. However, for rarer phenotypes, 
where neutral exploration is important, the T p are well 
approximated by Eq. (4). Again, the fraction of phenotypes 
that are immediately accessible should decrease for larger L. 
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5) In the intermediate regime /Z=10 , where NLfl=1.2, the 
population spreads over more phenotypes than in the 
monomorphic regime, but over fewer than in the polymorphic 
regime. Thus the crossover to the regime where neutral 
exploration is important occurs at a smaller (j) pq than for the 
monomorphic regime. In this intermediate fi regime neither 
Eq. (4) nor Eq. (7) suffices. Instead, we use the previously 
introduced factor y that multiplies N in Eq. (4) to achieve 
quantitative accuracy. In Appendix SI and in Figures S3 and 
S4, we explore the scaling of y with the parameters N,L,fi, for 
different dynamic regimes, and for a range of £. 

In summary then, our theory derived in the previous section 
accurately describes the median discovery time T p of this simple 
random GP map as a function of the parameters N,fl,tj) pq . We find 
that <j)„„!tiFp, and thus T p ~ 1 / F p in all regimes studied. The more 
frequent the phenotype, the earlier (and more often, see Figure S2) 
it appears as potentially selectable variation in an evolving 
population. Given the success of our theory for the random 
model, we now will test our theory and conclusions for a more 
complex GP map. 

RNA secondary structure mapping. One of the best 
studied GP mappings has RNA genotypes of length L made up 
of nucleotides G, C, U and A. The phenotypes are the minimum 
free-energy secondary structures for the sequences, which can be 
efficiently calculated [15]. The number of genotypes grows as 4 L , 
while the number of phenotypes is thought to grow roughly as 
Np~l.S L [4] so that Np«N(j. Moreover, sampling and exact 
enumerations [5, 16,22] have shown that the distribution of 
phenotype frequencies F p is highly biased, with a small fraction 
of phenotypes taking up the majority of genotypes. The neutral 
spaces Afq are typically broken up into a number of large 



components that are connected by single point mutations that 
allow neutral exploration [16,22]. By exhaustive enumeration of 
the L = 20 RNA mapping (see also Figure S5) we calculate the <j> pq 
between several neutral components of the 11,219 distinct 
secondary structures that the Ng = 4 20 x 1.1 x 10 12 genotypes 
map to. 

Figure 2b shows the $ pq for the largest component of the 
phenotype q drawn in the figure. This phenotype is ranked as the 
3rd most frequent for L = 20 and exhibits behaviour typical of this 
system. First, the <f> pq vary over many orders of magnitude. Second, 
as shown in the inset if <j> pq # 0, then the local </) pq are, to first order, 
proportional to the global F p . Finally, this neutral space connects to 
just over 75% of the total Np = 11,219 phenotypes in this 
particular map: Some <j> pq are zero even though F p can be quite 
large. Generally, the number of phenotypes that can be reached 
from Afq increases with F q [13,16]. 

We performed extensive simulations of the L = 20 RNA system. 
Typical results are shown in Figure 2B. First, we note that the 
median discovery times vary over many orders of magnitude. The 
most frequent are found in a median time of T p » 10 3 generations 
while after the maximum measured time of 2 x 10 9 generations, 
over 42% of the direcdy accessible phenotypes (with (j> pq # 0) have 
still not been found. We estimate that over 10 13 generations would 
be needed to discover all accessible phenotypes, giving a ten order 
of magnitude range in the T p . Second, the local frequency (f) pq is a 
good predictor for ranking T p (see Figure S6 for a comparison of 
T p and global frequency F p ). Further, the criterion (j) pq = 0 
accurately predicts that phenotypes are not discovered (see also 
Figure S6). However, in contrast to the random GP map, the T p 
are discovered at a slower rate than predicted by Eq. (7). Instead, 
we use a single y < 3Lp/N to renormalise N in Eq. (4). This slower 




Local frequency log 10 ^ Local phenotype rank Local frequency log 10 ^ 



Figure 2. Test of the meanfield model. A) Median discovery times T p for the random GP map averaged over 100 simulations with 7V= 1000 and 
varying mutation rates. Note that the y-axis is scaled with /<. In the the polymorphic limit (/i = 10~ 2 ), Eq. (4) (dashed line) describes discovery times 
well for (j> n < l/(NLn). Phenotypes with larger (j> M are part of the standing variation typically found in the first generation (yellow dash-dotted line). 
In the monomorphic limit (n = 10~ 6 ), Eq. (7) (dotted line) quantitatively describes T p for (l> M «(j>i, whereas Eq. (4) tracks the simulation data with just 
one fit parameter y = 0.099 multiplying N for the intermediate regime with n= 10~ 4 (solid line). For (/> s (j> L the curves follow Eq. (4), for reasons 
described in the text. Inset: For the random GP map the local phenotype frequency <j> pq correlates very well with the global frequency F p . B) Local 
frequency ij> pil ranked for the 8639 phenotypes that link with single point mutations from the \M q \ = 460,557,583 genotypes that map to this RNA 
structure; an example sequence from M q is shown in the figure. Inset: The local connections (j) pi] are roughly proportional to the global frequency F p , 
but there is significant scatter due to the internal correlations of the RNA neutral spaces. Organge points depict the 2580 phenotypes for which 
^ M =0. Light blue points depict the 4933 phenotypes that are discovered in our simulations, and the dark blue points depict the 3705 accessible 
phenotypes that are not found (q itself is shown in green). Q Simulations of T p (blue dots) versus (j) piJ for the RNA phenotype shown in B), compared 
to Eq. (4) (solid line) with a factor y = 0.070 multiplying N. Here 7V= 100, /< = 10~ 5 and the simulations were run for 2 x 10 9 generations. Also shown 
are the purely polymorphic (dashed) and monomorphic (dotted) predictions. Dark blue dots above 2 x 10 9 (dot-dashed line) depict some of the 3705 
accessible phenotypes that are not found (as can be seen in see the inset of B). We estimate that about 10° generations would be needed to find the 
phenotypes with the smallest ij> pq ^0. 
doi:1 0.1 371 /journal.pone.0086635.g002 
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discovery rate reflects the internal structure of the RNA: similar 
genotypes typically have similar mutational neighbourhoods [23], 
and so the population needs to neutrally explore longer in order to 
find novelty. Nevertheless, a single y factor yields a remarkably 
good fit for all the different phenotypes p (something we find for all 
source phenotypes q we have so far studied). Finally, we note that 
the three most frequent phenotypes are found relatively faster 
because they satisfy (j) pq s <j> L . As expected, for this larger system 
the fraction of phenotypes for which this holds is lower than for the 
random GP map with smaller L. 

Overall, the evolutionary dynamics of this rather complex RNA 
system resembles that of the much simpler random GP map. Most 
importandy, the discovery times vary over many orders of 
magnitude. More precisely, as long as <^, ? #0, T p ccl/(j> pq for 
both the monomorphic and polymorphic regimes: Phenotypic bias 
leads to a simple, systematic ordering in the discovery of novel 
phenotypes. 

The arrival of the frequent 

The many orders of magnitude difference in the arrival rate of 
variation between phenotypes should have many important 
implications for evolutionary dynamics. Consider for example 
the situation where the population has equilibrated to a phenotype 
q, which was the fitness peak, when subsequently the environment 
changes so that a different phenotype p has a higher fitness 1 + s. 
In order to fix, the alternative phenotype must first be found. If the 
time-scale Te on which the environment changes again is much 
longer than T p then it likely that the population will discover and 
fix p. However, if Te« T p , then a new phenotype p' may become 
more fit before p has time to fix. T p can vary over many orders of 
magnitude, so many potentially highly adaptive phenotypes may 
satisfy T p > Te and thus never be found. 

Consider also the situation where two phenotypes p\ and p2 are 
both more fit than q after an environmental change. If 
S2>si s 1/2N, then in a standard population genetics picture, 
we would expect pi to fix rather than pi as long as T p 2 < Te. 
However, this argument ignores the rate at which variation arises. 
If, for example, tj> p \ q » tj> p 2 q , then pi may fix well before p2 is 
discovered and fixes. 

To illustrate this effect, we study the L=12 RNA system 
depicted in Figure 3, where the source neutral space has 
A 9 = 1932 genotypes, while the two target phenotypes have 
<^, l9 = 0.067 and ^2^ = 0-0015, so (j> p2q /(j> plq x0.022, a relatively 
modest ratio compared to the what could be found from e.g. Fig. 2. 
For this particular system (j> p i p 2 = 0: there are no direct single 
mutation connections between the two target phenotypes - pi and 
p2 are distinct peaks of the fitness landscape. 

We simulated a population of N= 1000 individuals with fixed 
S\ =0.002> 1/2N, but with varying ratios Si/si > 1. The popu- 
lation begins on phenotype q and evolves until pi or p2 fixes. 

Results are shown in Figure 4. As the mutation rate increases 
and the system moves from the monomorphic to polymorphic 
regime, the probability that p2 is discovered at least once increases 
(and is largely independent of fitness). Nevertheless, phenotype pi 
is discovered much earlier and also much more often because 
<t> p\q y> § plq- Furthermore, in the monomorphic regime where £,» 1 
the population remains on a single genotype g much longer than it 
takes to explore all the neighbours. Thus if pi is accessible from g, 
then pi is likely arise repeatedly in relatively quick succession (in 
"bursts"). This effect, which arises naturally in our microscopic 
model [24], can significandy enhance the probability of fixation 
over that predicted by origin-fixation models [25] which ignore the 
discreteness of the source neutral space. 



Overall, our simulations show how the more frequent pheno- 
type pi can fix at the expense of the more fit phenotype p2. Given 
the many orders of magnitude difference possible between the T p , 
such an "arrival of the frequent" effect may prevent the arrival of 
the fittest: If a highly beneficial phenotype is never discovered, a 
much less adaptive but easily accessible phenotype may go to 
fixation instead. 

Finally, phenotype p2 is significantly less mutationally robust 
than pi (more frequent phenotypes are typically more robust 
[13,16]), and so once discovered, produces deleterious mutants at 
a higher rate, making it harder for p2 to fix at higher mutations 
rates, a phenomenon known as "survival of the flattest" [26], 
observed here for the lower ratios Sj/si at higher p.. Thus both the 
"arrival of the frequent" and the "survival of the flattest" mitigate 
against the fixation of phenotypes with lower frequency F p , even if 
their fitness is much higher. 

We note that differences in neutral network size have 
traditionally also been taken into account in terms of free fitness 
[27], which - in analogy with free energy in statistical physics [28] 
- incorporates an entropy-like component to account for 
mutational effects such as genetic drift and mutational robustness. 
This picture provides a theoretical foundation for the "survival of 
the flattest" [26] effect we observe at high mutation rates in 
Figure 4. However, the "arrival of the frequent" effect is 
fundamentally different because it does not rely on mutation- 
selection balance and quasi-equilibrium or steady-state assump- 
tions like free-fitness theory does. Rather, it reflects the strongly 
non-equilibrium effect that P2 is rarely or never found. In the 
example above, the difference in discovery times between p\ and 
P2 is rather modest, and so at large enough mutation rates P2 is 
found fairly regularly and free-fitness could be used to analyse 
results in that regime. But as can be seen for instance in Figure 2 
for L = 20 RNA, differences in discovery times can vary over many 
more orders of magnitude than is the case for our particular 
example, so that in practice highly adaptive yet rare phenotypes 
may not be discovered at all, even on very long timescales. 

Discussion 

Mutations provide the fuel for natural selection. Based on this 
principle, we have presented a detailed model of evolutionary 
dynamics that focuses on a microscopic description of the outcome 
of mutations. The phenotypic effect of mutations is mediated by 
the genotype-phenotype (GP) map which is therefore a crucial 
ingredient. As oudined in the introduction, several generic features 
are shared by many different example maps, independent of 
model details. Here we mainly focussed on the fact that these 
mapping are highly biased: Some phenotypes are realised by orders 
of magnitude more genotypes than most other phenotypes. 

Our calculations for a simplified random mapping and for the 
more complex RNA secondary structure model predict that the 
large bias observed in the GP maps translates into a similar order 
of magnitude variation in the median discovery times T p for a 
range of population genetic parameters. For both maps the local 
frequencies <j) pq (which predict discovery times) are a good 
predictor for the discovery times T p . For the random GP map 
<j> pq xF p . For RNA this relationship provides a rough first order 
estimate, but the local frequencies can also deviate strongly, 
especially when (j) pq = 0, which can occur even when the global 
frequency F p is large. For both maps the strong bias in the GP 
map leads to a systematic ordering of the median discovery times of 
alternative phenotypes, an effect that we postulate may hold for 
other GP maps as well. 
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Relative fitness 




Figure 3. Interconnections of neutral spaces in RNA influence evolutionary trajectories. A) £ = 12 RNA neutral component for phenotype 
q with jV (/ = 1932 genotypes (drawn in blue). Lines depict single mutations to itself, or to two alternative phenotypes pi (grey) and p2 (red). The 
genotypes were ordered using the Fruchterman-Reingold algorithm [30]. 6) Illustration of the fitness landscape. 
doi:1 0.1 371 /journal.pone.0086635.g003 



In light of the simplicity of our mean-field approximation, its 
success in predicting the first-discovery time T p (cf. Figure 2) is 
rather striking. In the random GP map, the excellent agreement 
probably arises because all genotypes in the source neutral space 
are similar in the sense that they have the same probability 
distribution to have a certain mutational neighbourhood. There 
are static fluctuations because the number of neighbours is less 
than the number of states with <j> ^ 0. But while these fluctuations 
have an effect on processes like fixation, they average out over the 
many runs used to find the mean or median T„ . By contrast, in the 
RNA GP map mutational neighbourhoods of adjacent genotypes 
are often correlated [13,23] so that a single neutral mutation does 
not completely re-shuffle the accessible phenotypes (as the mean- 
field assumption would assume). This effect explains why the value 
of the exploration parameter y we obtain by fitting is below the 
value suggested by our mean-field model, and also why we still 
observe around 1 order of magnitude variation in T p for very 
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Figure 4. The arrival of the frequent. Probability that phenotype p2 
is discovered (dotted lines) or is fixed (dashed lines) as a function of 
mutation rate p for different relative selection coefficients s 2 /s\ for 
Ns\ =2. The probability that p2 is discovered is independent of relative 
fitness (within statistical simulation errors). Phenotype pi is much more 
likely to fix than phenotype p2, even when the latter is much more fit, 
due to an "arrival of the frequent" phenomenon. 
doi:10.1371/journal.pone.0086635.g004 



similar values of (f> pq (see Figure 2). Despite such correlations 
(which we postulate may occur in other realistic GP maps), rare 
phenotypes (low <j) pq ) remain hard to find; the strong phenotypic 
bias in the RNA GP map provides a good a posteriori justification 
for our mean-field calculations: The many orders of magnitude 
range in (j) pq dominates the scale of the phenotype discovery times. 

The large differences we observe in the rate with which potential 
variation appears should have many consequences for evolutionary 
dynamics. There is of course a long history of invoking processes that 
impose directionality on the pathways available for evolutionary 
exploration (see ref. [29] for a recent discussion). Here, by solving 
microscopic population genetic models, we show in detail just how 
strong these orienting processes can be. Other authors have also 
pointed out how evolution may favour phenotypes with large neutral 
networks for RNA, see e.g. refs. [5,22]. Similar points have been 
made for protein models [12]. Consider, for example, our £ = 20 
RNA system. Despite its rather modest size, we find 10 orders of 
magnitude difference between the discovery times of frequent and 
rare phenotypes. These differences should be even more pronounced 
for larger £. In nature, selectable RNA phenotypes are of course 
characterised by more than just their secondary structure, and 
evolutionary processes don't always work at constant £. Neverthe- 
less, it is hard to see how such enormous variations in T p would not 
persist in some form in much more sophisticated treatments of 
biological RNA. Similar arguments can be made for the other GP 
maps we listed above. More generally we emphasise that including 
the GP map in population genetic calculations may be of importance 
to a wide range of evolutionary questions. 

We explicitly showed how phenotypes with a high local 
frequency can fix at the expense of locally rare phenotypes, even 
if the latter have much higher fitness. Taken together, these 
arguments suggest that the vast majority of possible phenotypes 
may never be found, and thus never fix, even though they may 
globally be the most fit: Evolutionary search is deeply non-ergodic. 
When Hugo de Vries was advocating for the importance of 
mutations in evolution, he famously said "Natural selection may 
explain the survival of the fittest, but it cannot explain the arrival 
of the fittest" [2] . Here we argue that the fittest may never arrive. 
Instead evolutionary dynamics can be dominated by the "arrival 
of the frequent". 
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Methods 

Simulations 

In the dynamic simulations, all N individuals of the population 
are initially assigned to a single random genotype in the source 
neutral space. Then the population evolves for 10JV generations to 
reach a steady-state dispersal on the neutral space before 
measurements are started. 

RNA 

Secondary structures for RNA were predicted from sequence 
using the Vienna package [15], version 1.8.5 with all parameters 
set to their default values. 

Supporting Information 

Appendix SI 

(PDF) 

Figure SI Static properties of the random GP map. A) Global 
phenotypes frequencies. In addition to the distribution of frequen- 
cies F p used in our simulations (orange), the diagram also shows the 
frequencies of RNA secondary structures at L=12, obtained by 
exhaustive enumeration using the Vienna package, Version 1.8.5 
with all parameters set to their default values [15]. B) Comparison of 
global frequencies F p and local frequencies <j> pq for the source 
neutral space q with rank 1 . The robustness of phenotype q (p = (j) qq ) 
is marked in green; alternative phenotypes (p¥^q) are shown in light 
blue. The dashed line marks the equality of global and local 
frequency F p =(j> pq . The relative size of deviations becomes more 
severe as F p becomes small: The less genotypes map into p, the less 
will frozen fluctuations in the GP map average out. 
(TIF) 

Figure S2 Total number of mutants per phenotype in 
different dynamic settings. The diagram shows the total 
number of mutants M p = Ym= \ m />(0 carrying phenotype p that 
were produced during a total of T=W 4 /(Np.) generations of 
simulation under the random GP map. Dots show the average 
over 100 simulations, error bars show the standard deviation. The 
dashed lines correspond to the mean-field theory M p = NLp(j> pq T 
that follows directiy from Eq. (3). In panels B and D, the 
populations are in the highly polymorphic regime (NLfi»l) and 
hence evolve towards greater robustness [19] so that the total 
number of non-neutral mutants is reduced. 
(TIF) 

Figure S3 Scaling of Ny with population dynamic 
parameters. The diagram shows the dependence of y on: A) 
mutation rate fi, B) population size N and £) number of mutants 
per generation NLfi. Note that the y-axis has been scaled by 
population size N. 
(TIF) 

Figure S4 Scaling of y with population dynamic param- 
eters. The diagram shows the dependence of y on: A) mutation 
rate fi, B) population size N and C] number of mutants per 
generation NLp,. In contrast to Figure S3, the y-axis shows y 
without any scaling factors. 
(TIF) 
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Figure S5 Phenotypic bias for RNA secondary struc- 
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tracing out all possible neutral mutations and counting how often 
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Figure S6 Predictions based on global frequency. The 

diagram shows the same median discovery times of alternative 
RNA secondary structures that are displayed in Figure 2c, but 
here as a function of the phenotypes' global frequencies F p rather 
than their local frequencies (j> pq . The different colors indicate: 
Accessible phenotypes that are typically discovered within the 
simulation time (T p (l/2)<2 x 10 9 ), <j> pq >0, light blue); accessible 
phenotypes that are typically not discovered (7^,(1/2) >2x 10 , 
(f> pq >0, dark blue); inaccessible phenotypes that are typically 
discovered {T p (\/2)<2 x 10 9 , <j> pq = 0, orange); inaccessible phe- 
notypes that are typically not discovered (7^,(1/2) > 2 x 10 , 
(j> pq = 0, red). The lines correspond to the prediction for T p based 
on global rather than local frequencies: T p (\/2) = \og2/(NLj.iF p ) 
(cf. Eq. (4)), dashed) and r / ,(l/2) = log2/(3L 2 W />) (cf Eq. (7)). 
In contrast to the predictions based on the local frequencies (j> pq in 
Figure 2c, we note the following: 1) Several phenotypes arise even 
earlier than predicted by the analogue of the polymorphic limit 
(points below dashed line). 2) Many phenotypes are not discovered 
even though other phenotypes of comparable (and even much 
lower) frequency do arise during the simulation. 3) 4 of the most 
frequent, but locally inaccessible phenotypes are discovered on a 
time-scale when double mutations become relevant (orange dots; 
since A =100 and fi= 10~ 5 , double mutants occur on the 
timescale t% ~ \/(N(Lp) 2 = 2.5 x 10 5 , so if double mutations were 
to lead to globally random phenotypes, we expect phenotypes with 
<l> pq = 0 to be discovered around T p xt2\og2/F p .) 
(TIF) 
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